/*
LP Solver based on GLPK Kit

*/


#include <glpk.h>

#ifndef PREREQ_H
#define PREREQ_H 
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <cmath>
#include "Timer.h"


#define index unsigned int
#define integer int
#define real float
#define dreal double
#define matrix_type float


#define INFTY 1.7e308

#define VERBOSE 1






/* ************************************************************************************************  *
                                  Reading and writing data
*  ************************************************************************************************* */

template<class T>
void printVector(T  *ptr, index size){
/*
	template prints array elements
*/
     index i;
     for(i=0;i<size;i++) std::cout  << ptr[i] << "  " ;
     std::cout << std::endl;
}

template<class T> 
void printMatrix(T *ptr, index rows,index cols)
{
/*
	template prints out matrix elements
*/
   for(index i=0;i<rows;i++)
   {
	   for(index j=0;j<cols;j++) std::cout << std::setw(3) << ptr[i+j*rows];
       std::cout << std::endl;
   }
     
   
}
template<class T>
long countNZ(T *ptr, size_t size)
{
	long nz = 0;
	for(index i=0;i<size;i++)
			nz += ptr[i]!=0 ? 1 : 0;
	
	
	return nz;
}




long getFileInput(const char *fileName,matrix_type** A, real **b, real **c, int &m, int &n){
/*
	prerequsities: A must correspond to the packing LP matrix
	reading from the file:
		1. read the dimension: m rows, n columns, number of non-zero elements 
		2. read matrix A as an array in column-wise fashion
		3. procedure reads number of non-zero entries for matrix A, then n elements for c and m elements for b
	
*/   
    FILE *filePt;      
    filePt = fopen(fileName,"r");
	long nonZeros;
    
   if(!filePt)
   { 
      printf("Input file missing ...\n"); 
	  exit(0);
   }          
   fscanf(filePt,"%d %d %ld",&m,&n,&nonZeros);    
       
   // memory allocation for arrays 
   *A = (matrix_type *) calloc(m*n,sizeof(matrix_type)); // matrix A
   *b = (real *) calloc(m,sizeof(real));
   *c = (real *) calloc(n,sizeof(real));
   
   index row = 0,col = 0;   
   integer readings = 0;      
   
   // read matrix A
   while(readings++<nonZeros)
   {
		matrix_type value; 
		
		fscanf(filePt,"%d %d %f",&row,&col,&value);
		(*A)[row+col*m] = (matrix_type) value; // column-wise storing
		
   }    
   // read array c
   readings = 0;
   while(readings<n)
   {
		real value;
		fscanf(filePt,"%f",&value);
		(*c)[readings] = (real) value;
                readings++;
   }
   // read array b
   readings=0;
   while(readings<m)
   {
		real value;
		fscanf(filePt,"%f",&value);
		(*b)[readings] = (real) value; 
		readings++;
   }
       
    fclose(filePt);
	
	return nonZeros;
  
};

#endif

int GLPKSol(matrix_type *A,real *b, real *c, dreal *x, dreal *y, integer m, integer n, dreal omega)
{
	if(VERBOSE)
	{
		std::cout << std::endl << std::endl;
		std::cout << "--------------------------------------------------------------------------" << std::endl;
		std::cout << "                             GLPK LP Solver                               " << std::endl;
		std::cout << "--------------------------------------------------------------------------" << std::endl;	
	}

	dreal primal, dual;

	long nz = countNZ(A,m*n);
	
	Timer compute_time;

	compute_time.start();
	glp_prob *lp_opt;
    glp_prob *lp_approx;
	long iter;


	if(VERBOSE) std::cout << "\n\rPreprocessing step ... " << std::endl;
	int *ia = new int[nz+1];
	int *ja = new int[nz+1];

	double *ar = new double[nz+1];
	double z_opt;
  
	ia[0] = 0; ja[0] = 0; ar[0]=0;
  
	index k=1;

	for (index i=0;i<m; i++)
		for(index j=0;j<n; j++)
	{
		if(abs(A[i+j*m])>0)
		{
			ia[k] = i+1;
			ja[k] = j+1;
			ar[k] =  abs(A[i+j*m]);
			k++;
		}
		else continue;
	}
	std::cout << "... done!" << std::endl;
		
	free(A);
  
    
  /* 
  define two LPs: optimal value LP, approx. value LP
  find optimal solution -> lp_opt
  */
  
  lp_opt = glp_create_prob();
  lp_approx = glp_create_prob();
  
  // params
  glp_smcp params1;
  glp_init_smcp(&params1);
  params1.msg_lev = GLP_MSG_ON;
  params1.presolve = GLP_ON;
  params1.meth = GLP_PRIMAL; 	

  // define constraints
  glp_set_obj_dir(lp_opt, GLP_MAX);  // set optimization problem as maximization
  glp_add_rows(lp_opt, m);           // define m rows
	
  for (int i=1; i<=m; i++) glp_set_row_bnds(lp_opt, i, GLP_UP, 0.0, b[i-1]); // set constraints
	
  glp_add_cols(lp_opt, n); // define number of columns
  for (int i=1; i<=n; i++)
  {
	glp_set_col_bnds(lp_opt, i, GLP_LO, 0.0, 0.0);  // unbounded for cols
	glp_set_obj_coef(lp_opt, i, c[i-1]); // define cost vector
  }
  
  // load constraint matrix into LP
  glp_load_matrix(lp_opt, nz, ia, ja, ar);

  // copy lp_opt to lp_approx
  glp_copy_prob(lp_approx,lp_opt,GLP_OFF);
  
  if(omega==0.0)
  {
	if(VERBOSE) std::cout << "Finding OPT value ..." << std::endl;
	 		  	
	
	glp_simplex(lp_opt, &params1);
	z_opt = glp_get_obj_val(lp_opt);
		
	primal = dual = z_opt;	
	// get the solution
	for (int i=1; i<= n; i++)
		x[i] = glp_get_col_prim(lp_opt, i);
	
	for (int i=1; i<= m; i++)
		y[i] = glp_get_col_dual(lp_opt, i);
		
	// get number of iterations
	 iter = lpx_get_int_parm(lp_opt, LPX_K_ITCNT);

	glp_delete_prob(lp_opt);
	
	}
	
	if (omega>0)
	{
		 
	    if(VERBOSE) std::cout << "Finding aapproximate solution  to LP ... " << std::endl;
		std::cout << "Procedure not yet implemented ...";

		glp_delete_prob(lp_approx);
		return 1;      			
		
	}

	
	
	compute_time.stop();
	
	if(VERBOSE)
	{
    std::cout << "\nDATA SET: rows: " << m << ", cols: " << m << ", density: " << nz/(double)(n*m)*100<< "%" << std::endl;
	std::cout << "---------------------------------------------------------------" << std::endl;
			
	std::cout << "Number of iterations:                           | " << iter << std::endl;
	std::cout << "----------------------------------------------------------------" << std::endl;
	std::cout << "*** Compute time:                               | " << compute_time.secs() << std::endl;
	std::cout << "----------------------------------------------------------------" << std::endl;
	
	printf(" primal value: %f       \n\r dual value: %f     \n\r approx. ratio: %f \n",primal,dual,dual/primal);
	std::cout << "---------------------------------------------------------------------------------------" << std::endl;
	}

		
	delete [] x;
	delete [] y;
	
	
	return 1;

}